source("/users/brucedesmarais/Dropbox/SenPE/Code/CreateSenPENets.R")

library(igraph)
community <- list()

for(t in 1:length(mats)){
mati <- mats[[t]]
ntype <- c(rep(0,nrow(mati)),rep(1,ncol(mati)))
el <- NULL
for(i in 1:nrow(mati)){
	for(j in 1:ncol(mati)){
		if(mati[i,j]==1){
			n1 <- i-1
			n2 <- nrow(mati)+j-1
			el <- rbind(el,c(n1,n2))
			}
		}
	}

# Project the bipartite network

g <- graph.bipartite(ntype,c(t(el)))
g2 <- bipartite.projection(g)[[2]]
m <- as.matrix(g2)
edge_g2 <- cbind(m[4][[1]],m[3][[1]])
wg2 <- numeric(nrow(edge_g2))
for(i in 1:nrow(edge_g2)){
		indi <- edge_g2[i,1]
		indj <- edge_g2[i,2]
		wg2[i] <- sum(mati[,indi+1]*mati[,indj+1])
	} 

cg2 <- fastgreedy.community(g2)
cg2w <- fastgreedy.community(g2,weights=wg2)
mem <- community.to.membership(g2,cg2$merges,which.max(cg2$modularity))
memw <- community.to.membership(g2,cg2w$merges,which.max(cg2w$modularity))
nett <- get.adjacency(g2)
nett[edge_g2+1] <- wg2
community[[t]] <- list(cg2,cg2w,mem,memw,nett)
}

fowler <- read.csv("SH.csv",stringsAsFactors=F)
fowler <- subset(fowler,fowler[,2]=="S")


comid <- numeric(nrow(id_master))

for(t in 4:length(mats)){
	congi <- congs[t]
	ind <- match(paste(colnames(mats[[t]]),congi,sep=""),paste(id_master$sen_pe,id_master$Congress,sep=""))
	comid[ind] <- community[[t]][[4]]$membership
	}
	
comid[which(id_master$Congress>105)] <- NA
	
write.csv(data.frame(id_master,comid),file="community_membership.csv")

## Get Modularity By Both Community and Party ##

mod_mat <- NULL

for(t in 4:length(mats)){
mati <- mats[[t]]
ntype <- c(rep(0,nrow(mati)),rep(1,ncol(mati)))
el <- NULL
for(i in 1:nrow(mati)){
	for(j in 1:ncol(mati)){
		if(mati[i,j]==1){
			n1 <- i-1
			n2 <- nrow(mati)+j-1
			el <- rbind(el,c(n1,n2))
			}
		}
	}
	

# Project the bipartite network

g <- graph.bipartite(ntype,c(t(el)))
g2 <- bipartite.projection(g)[[2]]
m <- as.matrix(g2)
edge_g2 <- cbind(m[4][[1]],m[3][[1]])
wg2 <- numeric(nrow(edge_g2))
for(i in 1:nrow(edge_g2)){
		indi <- edge_g2[i,1]
		indj <- edge_g2[i,2]
		wg2[i] <- sum(mati[,indi+1]*mati[,indj+1])
	} 

mem <- community[[t]][[4]]$membership
maxMod <- modularity(g2,mem,wg2)
memsi <- colnames(mats[[t]])
congi <- congs[t]
idmi <- subset(id_master,id_master$Congress==congi)
fowi <- subset(fowler,fowler$congress==congi)
fowid <- idmi$cosp_summary[match(memsi,idmi$sen_pe)]
party <- fowi$party[match(fowid,fowi$labels)]
party <- as.numeric(party==200)
party[which(is.na(party))] <- 0
partyMod <- modularity(g2,party,wg2)

mod_mat <- rbind(mod_mat,c(congi,maxMod,partyMod))

}

## Plot the degree distribution, with mean and median
detach(package:igraph)
library(sna)
for(i in 2:length(wnets)){
filei <- paste("/users/brucedesmarais/Dropbox/SenPE/Tex/degdist",names(wnets)[i],".pdf",sep="")
pdf(filei,height=4.5, width=5, family="Times", pointsize=13)
par(las=1,mar=c(6,6,.5,.5),cex.lab=2,cex.axis=1.75)
wni <- wnets[[i]] > 0
degi <- degree(wni,gmode="graph")
hist(degi,col="grey70",border="grey45",main="",xlab="",xlim=c(0,80),ylim=c(0,40),breaks=c(0:8)*10)
abline(v=c(mean(degi),median(degi)),lwd=2.75,lty=c(1,2))
title(xlab="Degree",line=3.5)
dev.off()
}
#Fowler's Data, Dem = 100, Rep = 200

fowler <- read.csv("SH.csv",stringsAsFactors=F)
fowler <- subset(fowler,fowler[,2]=="S")

# Plot communities based on party composition
library(colorRamps)
pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/communities.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(96:105,seq(0,100,length=length(4:13)),col="white",xlab="",ylab="Cumulative Percent",bty="n",xlim=c(95.5,105.5),xaxt="n")
red2blues <- gray((100:1)/100)

for(i in 4:13){
	memsi <- colnames(mats[[i]])
	csizei <- community[[i]][[4]]$csize
	cmemi <- community[[i]][[4]]$membership
	coms <- 0:(max(which(csizei>1)-1))
	congi <- congs[i]
	csizei <- csizei[which(csizei>1)]
	props <- 100*csizei/sum(csizei)
	pdem <- numeric(length(coms))
	fowi <- subset(fowler,fowler$congress==congi)
	idmi <- subset(id_master,id_master$Congress==congi)
	for(j in 1:length(pdem)){
		cmemj <- memsi[which(cmemi==coms[j])]
		fowid <- idmi$cosp_summary[match(cmemj,idmi$sen_pe)]
		party <- fowi$party[match(fowid,fowi$labels)]
		pdem[j] <- mean(party==100,na.rm=T)
		}
	cszsrt <- csizei[order(pdem)]
	pdem <- pdem[order(pdem)]
	cols <- red2blues[1+round(pdem*99)]
	cssz <- cumsum(cszsrt/sum(cszsrt))
	lb <- c(0,cssz[1:(length(cssz)-1)])*100
	ub <- cssz*100
	for(j in 1:length(pdem)){
		rect(xleft=congi-.5,xright=congi+.5,ybottom=lb[j],ytop=ub[j],col=cols[j])
		}
	}

axis(1,at=96:105)
dev.off()

# Plot max modularity vs. polarization

pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/modularity.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(mod_mat[,c(1,2)],ylim=c(min(c(mod_mat[,c(2,3)])),max(c(mod_mat[,c(2,3)]))), lwd=2.5,xlab="",ylab="Modularity",type="l")
lines(mod_mat[,c(1,3)],col="grey55",lwd=2.5)
abline(v=96:105,lwd=2.5,col="grey70",lty=3)
points(mod_mat[,c(1,2)],cex=1.5,pch=19)
points(mod_mat[,c(1,3)],pch=23,cex=1.5,col="grey55",bg="grey60")
legend("bottomright",legend=c("Max Modularity","Polarization"),col=c("black","grey55"),lwd=2.5,pch=c(19,23),pt.bg=c("black","grey55"),bg="white")

axis(1,at=96:105)
dev.off()

# Clustering Distance

dclust <- function(mem,net){
	umem <- unique(mem)
	dmat <- matrix(0,length(umem),length(umem))
	dnet <- matrix(0,nrow(net),nrow(net))
	for(i in 1:length(umem)){
		for(j in 1:length(umem)){
			if(i != j){
				ci <- which(mem==umem[i])
				cj <- which(mem==umem[j])
				dmat[i,j] <- 1/(.1+mean(net[ci,cj]))
				dnet[ci,cj] <- 1/(.1+mean(net[ci,cj]))
				}
			}
		}
	dnet
	}

corLab <- read.csv("/Users/brucedesmarais/Dropbox/SenPE/Data/PartyCorrect.csv",stringsAsFactors=F)

cinds <- match(97:105,congs)
for(c in cinds){
filec <- paste("/users/brucedesmarais/Dropbox/SenPE/Tex/CCirc",congs[c],".pdf",sep="")
pdf(filec,height=6, width=6, family="Times", pointsize=8.5)
par(las=1,mar=c(2,2,2,2),xpd=T)

mems <- community[[c]][[4]]$membership
net <- community[[c]][[5]]
clustd <- dclust(mems,net)
colnames(clustd) <- colnames(mats[[c]])
rownames(clustd) <- colnames(mats[[c]])

labs <- colnames(clustd)
if(c==11){
for(i in 1:length(labs)){
	labi <- labs[i]
	vlab <- strsplit(labi,split="")[[1]]
	com <- which(vlab==",")
	if(length(com) >0){
	if(nchar(labi) > com+2){
	labi <- substr(labi,1,com+2)
	}
	}
	labs[i] <- labi
	}
}

colnames(clustd) <- labs
rownames(clustd) <- labs

hc <- hclust(as.dist(clustd))

labs <- colnames(clustd)
if(c==11){
for(i in 1:length(labs)){
	labi <- labs[i]
	vlab <- strsplit(labi,split="")[[1]]
	com <- which(vlab==",")
	if(length(com) >0){
	if(nchar(labi) > com+2){
	labi <- substr(labi,1,com+2)
	}
	}
	labs[i] <- labi
	}
}

#sen_dat <- read.dta("fowlersenpenbinomdata.dta")
sen_dat <- read.dta("senator_data.dta")


cols <- rep("black",nrow(clustd))
ind <- match(labs,sen_dat$name)
isna <- which(is.na(ind))
ind[isna] <- 91
party <- sen_dat$party[ind]
party[isna] <- corLab$Party[match(labs[isna],corLab$Name)]
cols[which(party=="R")] <- "red"
cols[which(party=="D")] <- "blue"


# Circular dendrogram in R	 	 
library(ape)
library(cluster) 
data(mtcars)
plot(as.phylo(hc),type="fan",cex=1.25,edge.width=2,font=1,label.offset=.25,tip.color=cols)	 
text(0,-1.5,lab=congs[c],cex=3)
dev.off()
}


